#ifndef OKADA_HXX
#define OKADA_HXX

#include "FTensor.hpp"

class Okada
{
public:
  static const double eps;
  const double alpha, width, length,
    alp1, alp2, alp3, alp4, alp5;
  double cos_dip2, sin_dip2, sin_cos_dip;
  FTensor::Tensor1<double,3> dislocation, origin;
  FTensor::Tensor2_symmetric<double,3> rotation_strike, z_transform_dip;
  FTensor::Tensor2<double,3,3> rotation_dip;

  double radians(const double &degrees)
  {
    const double pi(4*atan(1.0));
    return degrees*pi/180;
  }

  double clamp(const double &value)
  {
    if(std::abs(value)<eps)
      return 0;
    return value;
  }

  Okada(const double &lambda, const double &mu, const double &slip,
        const double &Opening,
        const double &Width,
        const double &Length,
        const double &Strike,
        const double &Dip,
        const double &Rake,
        const FTensor::Tensor1<double,3> &Origin);

  typedef std::pair<FTensor::Tensor1<double,3>,FTensor::Tensor2<double,3,3> >
  Displacement;
  Displacement displacement(const FTensor::Tensor1<double,3> &coord);

  Displacement dc3d(const FTensor::Tensor1<double,3> &coord);

  Displacement source_contribution(const bool &real_source, const double xi[],
                                   const FTensor::Tensor1<double,3> &coord);
};

#endif
